/*
 * case.cpp
 *
 *  Created on: 11-09-2012
 *      Author: mith
 */

#include "case_BF_sW.h"
// MY HEADERS
#include"mros-v42.h"

namespace case_BF_sW
{
debug deb;
global glob;

// *******************************************************************************************  //
// *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY ***
// *******************************************************************************************  //

#define Sqrt sqrt
#define Power pow

/*
 *********************************************************
 ******************** Structures *************************
 *********************************************************
*/

struct parameters {
	double n;
	double t;
	double J;
	double U;
	double beta;
	double alpha;
	unsigned int L;
};

struct variables {
	double chiS;
	double DeltaS;
	double chiB;
	double DeltaB;
	double mu;
	double lBF;

	// extra
	double Energy;
	double nB;
 };

/*
 *********************************************************
 ***************** Small useful functions ****************
 *********************************************************
*/

inline void clear(parameters &p) {
	p.n = p.t = p.J = p.U = p.beta = p.alpha = 0.;
	p.L = 0;
}

inline void clear(variables &v) {
	v.chiS = v.DeltaS = v.chiB = v.DeltaB = v.mu = v.lBF = 0.;
	v.Energy = v.nB;
}

inline void addTo(const variables &q, variables &v, double d=1.) {
	v.chiS += q.chiS*d;
	v.DeltaS += q.DeltaS*d;
	v.chiB += q.chiB*d;
	v.DeltaB += q.DeltaB*d;
	v.mu += q.mu*d;
	v.lBF += q.lBF*d;
	v.Energy += q.Energy*d;
	v.nB += q.nB*d;
}

inline void copyTo(variables &v, const variables &q, double d=1.) {
	v.chiS = q.chiS*d;
	v.DeltaS = q.DeltaS*d;
	v.chiB = q.chiB*d;
	v.DeltaB = q.DeltaB*d;
	v.mu = q.mu*d;
	v.lBF = q.lBF*d;
	v.Energy = q.Energy*d;
	v.nB = q.nB*d;
}

inline void copyTo(parameters &v, const parameters &q) {
	v.n = q.n;
	v.t = q.t;
	v.J = q.J;
	v.U = q.U;
	v.beta = q.beta;
	v.alpha = q.alpha;
	v.L = q.L;
}

inline double abs(double a) {
	if(a<0) a=-a;
	return a;
}

/*
 *********************************************************
 ***************** CLASS Functions ***********************
 *********************************************************
*/

class Func {
private:
	parameters p; // parameters
	variables v; // variables
	double chiS, DeltaS, chiB, DeltaB, mu, lBF, n; // contractions

	// Abbreviations - dependend on kx, ky
	double gammak;
	double etak;

public:
	double d, lchiS, lDeltaS, lchiB, lDeltaB, ln; // contractions

	double qt;
	double dqtdn;
	double dqtdd;
	double qs;
	double dqsdn;
	double dqsdd;

	double E1[2];
	double dE1dlchiS[2];
	double dE1dlDeltaS[2];
	double dE1dlchiB[2];
	double dE1dlDeltaB[2];
	double dE1dln[2];
	double dE1dlBF[2];

	double E2[2];
	double dE2dlchiS[2];
	double dE2dlDeltaS[2];
	double dE2dlchiB[2];
	double dE2dlDeltaB[2];
	double dE2dln[2];
	double dE2dlBF[2];

	Func(const parameters &p, const variables &v) {
		set(p, v);
		refresh(); // Było zakomentowane?
	}

	void set(const parameters &p, const variables &v) {
		copyTo(this->p, p);
		copyTo(this->v, v);
		refresh();
	}

	void set(const parameters &p) {
		copyTo(this->p, p);
		refresh();
	}

	// Gutzwiller
	double fqt();
	double fdqtdn();
	double fdqtdd();
	double fqs();
	double fdqsdn();
	double fdqsdd();

    // Energy
    double fE1(double gammak, double etak, int l);
    double fdE1dlchiS(double gammak, double etak, int l);
    double fdE1dlDeltaS(double gammak, double etak, int l);
    double fdE1dlchiB(double gammak, double etak, int l);
    double fdE1dlDeltaB(double gammak, double etak, int l);
    double fdE1dln(double gammak, double etak, int l);
    double fdE1dlBF(double gammak, double etak, int l);

    double fE2(double gammak, double etak, int l);
    double fdE2dlchiS(double gammak, double etak, int l);
    double fdE2dlDeltaS(double gammak, double etak, int l);
    double fdE2dlchiB(double gammak, double etak, int l);
    double fdE2dlDeltaB(double gammak, double etak, int l);
    double fdE2dln(double gammak, double etak, int l);
    double fdE2dlBF(double gammak, double etak, int l);


	// REFRESH
	void refresh() {
		//variable
		chiS = v.chiS;
		DeltaS = v.DeltaS;
		chiB = v.chiB;
		DeltaB = v.DeltaB;
		mu = v.mu;
		lBF = v.lBF;
		n = p.n;

		// obliczam funkcje Gutzwillara i poszczegolne pochodne
		qt		= fqt();
		dqtdn = fdqtdn();
		dqtdd = fdqtdd();
		qs		= fqs();
		dqsdn = fdqsdn();
		dqsdd = fdqsdd();

		// obliczam odpowiednie
		d = ( -8.*p.t*chiS*(chiB*chiB+DeltaB*DeltaB)*dqtdd + 3*(1-p.alpha)*p.J*(chiS*chiS+DeltaS*DeltaS)*dqsdd - p.alpha*p.J*(3.*chiB*chiB+DeltaB*DeltaB)*dqsdd )/(2*p.U);
		lchiS = -qt*p.t*(chiB*chiB+DeltaB*DeltaB) + 3.*(1-p.alpha)*qs*p.J*chiS/4.;
		lDeltaS = 3.*(1-p.alpha)*qs*p.J*DeltaS/4.;
		lchiB = -2.*qt*p.t*chiS*chiB - 3.*p.alpha*qs*p.J*chiB/4.;
		lDeltaB = -2.*qt*p.t*chiS*DeltaB - p.alpha*qs*p.J*DeltaB/4.;
		ln = -8.*p.t*chiS*(chiB*chiB+DeltaB*DeltaB)*dqtdn + 3*(1-p.alpha)*p.J*(chiS*chiS+DeltaS*DeltaS)*dqsdn - p.alpha*p.J*(3.*chiB*chiB+DeltaB*DeltaB)*dqsdn ;

		// zeby pozniej mozna bylo wypisac
		glob.d = d;
		glob.lchiS = lchiS;
		glob.lDeltaS = lDeltaS;
		glob.lchiB = lchiB;
		glob.lDeltaB = lDeltaB;
		glob.ln = ln;
		glob.qt = qt;
		glob.qs = qs;
		glob.Energy = v.Energy;		// potrzebne?
		glob.nB = v.nB;				// potrzebne?
	}

	//REFRESH XY
	void refresh(double kx, double ky) {
		gammak = 2.*(cos(kx) + cos(ky));
		etak = 2.*(cos(kx) - cos(ky));
		refresh();

		for(int i=0; i<2; i++) { // 2 families
			E1[i]			= fE1(gammak, etak, 1-2*i);
			dE1dlchiS[i]	= fdE1dlchiS(gammak, etak, 1-2*i);
			dE1dlDeltaS[i]	= fdE1dlDeltaS(gammak, etak, 1-2*i);
			dE1dlchiB[i]   	= fdE1dlchiB(gammak, etak, 1-2*i);
			dE1dlDeltaB[i]  = fdE1dlDeltaB(gammak, etak, 1-2*i);
			dE1dln[i] 		= fdE1dln(gammak, etak, 1-2*i);
			dE1dlBF[i] 		= fdE1dlBF(gammak, etak, 1-2*i);

			E2[i]			= fE2(gammak, etak, 1-2*i);
			dE2dlchiS[i]	= fdE2dlchiS(gammak, etak, 1-2*i);
			dE2dlDeltaS[i]	= fdE2dlDeltaS(gammak, etak, 1-2*i);
			dE2dlchiB[i]   	= fdE2dlchiB(gammak, etak, 1-2*i);
			dE2dlDeltaB[i]  = fdE2dlDeltaB(gammak, etak, 1-2*i);
			dE2dln[i] 		= fdE2dln(gammak, etak, 1-2*i);
			dE2dlBF[i] 		= fdE2dlBF(gammak, etak, 1-2*i);
		}
	}
};


/////////////////
/////////////////
/////////////////

double Func::fqt() {
	return (2*Power(Sqrt(Power(d,2)) + Sqrt(1 + Power(d,2) - n),2)*(2*Power(d,2) - n))/((-2 + n)*n);
}

double Func::fdqtdn() {
	return (-2*(Sqrt(Power(d,2)) + Sqrt(1 + Power(d,2) - n))*(-4*Power(Power(d,2),1.5)*Sqrt(1 + Power(d,2) - n) + 4*Power(d,4)*(-1 + n) + Power(d,2)*(-4 + (4 - 3*n)*n) +
		       n*(4*Power(Power(d,2),1.5)*Sqrt(1 + Power(d,2) - n) + n - Sqrt(Power(d,2))*Sqrt(1 + Power(d,2) - n)*n)))/(Sqrt(1 + Power(d,2) - n)*Power(-2 + n,2)*Power(n,2));
}

double Func::fdqtdd() {
	return (4*d*Power(Sqrt(Power(d,2)) + Sqrt(1 + Power(d,2) - n),2)*(2*Power(d,2) + 2*Sqrt(Power(d,2))*Sqrt(1 + Power(d,2) - n) - n))/(Sqrt(Power(d,2))*Sqrt(1 + Power(d,2) - n)*(-2 + n)*n);
}


double Func::fqs() {
	return (4*Power(-2*Power(d,2) + n,2))/(Power(-2 + n,2)*Power(n,2));
}

double Func::fdqsdn() {
	return (-8*(-2*Power(d,2) + n)*(-4*Power(d,2)*(-1 + n) + Power(n,2)))/(Power(-2 + n,3)*Power(n,3));
}

double Func::fdqsdd() {
	return (32*d*(2*Power(d,2) - n))/(Power(-2 + n,2)*Power(n,2));
}

////////////////////////////////////////////////////////////////////////////
double Func::fE1(double gammak, double etak, int l) {
	return -(l*Sqrt(Power(etak,2)*Power(lDeltaS,2) + Power(gammak*lchiS + lBF + ln + mu,2)));
}

double Func::fdE1dlchiS(double gammak, double etak, int l) {
	return -((gammak*l*(gammak*lchiS + lBF + ln + mu))/Sqrt(Power(etak,2)*Power(lDeltaS,2) + Power(gammak*lchiS + lBF + ln + mu,2)));
}

double Func::fdE1dlDeltaS(double gammak, double etak, int l) {
	return -((Power(etak,2)*l*lDeltaS)/Sqrt(Power(etak,2)*Power(lDeltaS,2) + Power(gammak*lchiS + lBF + ln + mu,2)));
}

double Func::fdE1dlchiB(double gammak, double etak, int l) {
	return 0;
}

double Func::fdE1dlDeltaB(double gammak, double etak, int l) {
	return 0;
}

double Func::fdE1dln(double gammak, double etak, int l) {
	return -((l*(gammak*lchiS + lBF + ln + mu))/Sqrt(Power(etak,2)*Power(lDeltaS,2) + Power(gammak*lchiS + lBF + ln + mu,2)));
}

double Func::fdE1dlBF(double gammak, double etak, int l) {
	return -((l*(gammak*lchiS + lBF + ln + mu))/Sqrt(Power(etak,2)*Power(lDeltaS,2) + Power(gammak*lchiS + lBF + ln + mu,2)));
}

////////////////////////////////////////////////////////////////////////////

double Func::fE2(double gammak, double eta, int l) {
	return -(gammak*(lchiB + l*lDeltaB)) + lBF;
}

double Func::fdE2dlchiS(double gammak, double etak, int l) {
	return 0;
}

double Func::fdE2dlDeltaS(double gammak, double etak, int l) {
	return 0;
}

double Func::fdE2dlchiB(double gammak, double etak, int l) {
	return -gammak;
}

double Func::fdE2dlDeltaB(double gammak, double etak, int l) {
	return -(gammak*l);
}

double Func::fdE2dln(double gammak, double etak, int l) {
	return 0;
}

double Func::fdE2dlBF(double gammak, double etak, int l) {
	return 1.;
}

/*
 *********************************************************
 ******************* CLASS Equations *********************
 *********************************************************
*/

class Equa {
private:
	parameters p; // parameters
	variables v; // variables
	Func *f;

	void onePoint(variables &y, int i, int j);

public:
	Equa(parameters &p, variables &v) : p(p), v(v) { f = new Func(p, v); }
	~Equa() { delete f; }

	void setPar(parameters &p) {
		this->p = p;
		f->set(p);
	}

	void setParVar(parameters &p, variables &v) {
		this->p = p;
		this->v = v;
		f->set(p, v);
	}

	void run2(variables &w);

};

void Equa::onePoint(variables &w, int i, int j) {
	double x = (double)i*2.*M_PI/p.L;
	double y = (double)j*2.*M_PI/p.L;

	f->refresh(x, y);
	//if(DEEP_DEBUG) { f->printLambda(); }

	double E1tan[2] = { (1-tanh(p.beta*f->E1[0]/2.))/2., (1-tanh(p.beta*f->E1[1]/2.))/2.};
	double E2tan[2] = { (1-tanh(p.beta*f->E2[0]/2.))/2., (1-tanh(p.beta*f->E2[1]/2.))/2.};

	for(int i=0; i<2; i++) {
		// up
		w.chiS		+= E1tan[i]*f->dE1dlchiS[i] 	+ E2tan[i]*f->dE2dlchiS[i];
		w.DeltaS	+= E1tan[i]*f->dE1dlDeltaS[i] 	+ E2tan[i]*f->dE2dlDeltaS[i];
		w.chiB 		+= E1tan[i]*f->dE1dlchiB[i] 	+ E2tan[i]*f->dE2dlchiB[i];
		w.DeltaB	+= E1tan[i]*f->dE1dlDeltaB[i] 	+ E2tan[i]*f->dE2dlDeltaB[i];
		w.mu 		+= E1tan[i]*f->dE1dln[i] 		+ E2tan[i]*f->dE2dln[i];
		w.lBF 		+= E1tan[i]*f->dE1dlBF[i] 		+ E2tan[i]*f->dE2dlBF[i];

		w.nB += E2tan[i];

        // it would be safer have following code - warning: 30 is chosen arbitrarily
		if(p.beta*f->E1[i] > 30.) { w.Energy += -exp(-p.beta*f->E1[i])/p.beta; } // basically 0
		else if(p.beta*f->E1[i] < -30) { w.Energy += f->E1[i]; }
		else { w.Energy += -log(1 + exp(-p.beta*f->E1[i]))/p.beta; }

		if(p.beta*f->E2[i] > 30.) { w.Energy += -exp(-p.beta*f->E2[i])/p.beta; } // basically 0
		else if(p.beta*f->E2[i] < -30) { w.Energy += f->E2[i]; }
		else { w.Energy += -log(1 + exp(-p.beta*f->E2[i]))/p.beta; }
	}
}

void Equa::run2(variables &w)
{
	clear(w);
	variables temp; clear(temp); // temporary vector using during summation
	int max = p.L/2; // attention: L should be even!

	//int iii = 0;

    // SERIAL CASE - OK
    for(int i=-max; i<max; i++) {
		for(int j=-max; j<max; j++) {
			onePoint(temp, i, j);
     	}
    }
    addTo(temp, w, 1);
    clear(temp);

    // SMART WAY - ALSO OK
    // SERIAL CASE
   	// SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
   	// * edges (4 times)
/*   	for(int i=1; i<max; i++) {
   		onePoint(temp, i, i);
   		onePoint(temp, max, i);
   		onePoint(temp, i, 0);
   		//iii += 12;
   	}
      addTo(temp, w, 4);
      clear(temp);
      // * interior area (8 times) - SERIAL CASE
      for(int i=1; i<max; i++) {
   		for(int j=1; j<i; j++) {
   			onePoint(temp, i, j); //iii+=8;									//	sleep(1);
        	}
   	}
   	addTo(temp, w, 8);
   	clear(temp);
   	// * missing corners
   	onePoint(w, max, max);
   	onePoint(w, max, 0);
   	onePoint(w, 0, max);
   	onePoint(w, 0, 0);
   	//iii += 4;
*/
    //printf("%d\n", iii);

	// what was out of sumation
	double L2 = (double)(p.L*p.L);
	w.chiS /= L2;
	w.DeltaS /= L2;
	w.chiB /= L2;
	w.DeltaB /= L2;
	w.mu /= L2;
	w.lBF /= L2;
    w.Energy /= L2;
    w.nB /= L2;

	w.chiS += 8.*v.chiS;
	w.DeltaS += 8.*v.DeltaS;
	w.chiB += 8.*v.chiB;
	w.DeltaB += 8.*v.DeltaB;
	w.mu += p.n - 1.;
	w.lBF += -2.;

    v.Energy = w.Energy + 8.*f->qt*p.t*v.chiS*(v.chiB*v.chiB+v.DeltaB*v.DeltaB) - 3.*(1-p.alpha)*f->qs*p.J*(v.chiS*v.chiS+v.DeltaS*v.DeltaS) +\
    		p.alpha*f->qs*p.J*(3.*v.chiB*v.chiB+v.DeltaB*v.DeltaB) + p.U*f->d*f->d + 8.*(f->lchiS*v.chiS + f->lDeltaS*v.DeltaS +\
    		f->lchiB*v.chiB + f->lDeltaB*v.DeltaB) + f->ln*(p.n-1.) - 2.*v.lBF - v.mu;
    v.nB = w.nB;

    // results
    case_BF_sW::glob.Energy = v.Energy + v.mu*p.n;
    case_BF_sW::glob.nB = v.nB;
}

int function(const gsl_vector * x, void *params, gsl_vector * f) {

	struct parameters p;
	p.n = ((Par *) params)->n();
	p.t = ((Par *) params)->t();
	p.J = ((Par *) params)->J();
	p.U = ((Par *) params)->U();
	p.beta = ((Par *) params)->beta();
	p.alpha = ((Par *) params)->alpha();
	p.L = ((Par *) params)->L();

	variables v;
	v.chiS = gsl_vector_get (x, 0);
	v.DeltaS = gsl_vector_get (x, 1);
	v.chiB = gsl_vector_get (x, 2);
	v.DeltaB = gsl_vector_get (x, 3);
	v.mu = gsl_vector_get (x, 4);
	v.lBF = gsl_vector_get (x, 5);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5));
	}

    // It's start of the calculations

    Equa E(p, v);
    variables w;
    //E.run(w);
    //E.run_global(w);
    E.run2(w);

    // It's end of the calculations

	gsl_vector_set (f, 0, w.chiS);
	gsl_vector_set (f, 1, w.DeltaS);
	gsl_vector_set (f, 2, w.chiB);
	gsl_vector_set (f, 3, w.DeltaB);
	gsl_vector_set (f, 4, w.mu);
	gsl_vector_set (f, 5, w.lBF);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5));
		printf("###  %24.15e# %25.15e# %24.15e# %24.15e# %24.15e# %24.15e# | E = %24.15e nB = %24.15e\n",
			gsl_vector_get (f, 0), gsl_vector_get (f, 1), gsl_vector_get (f, 2),
			gsl_vector_get (f, 3), gsl_vector_get (f, 4), gsl_vector_get (f, 5),
			case_BF_sW::glob.Energy, case_BF_sW::glob.nB);
		printf("\n");
	}

    return GSL_SUCCESS;
}

// ************* CLASSES ************* //

// Par

void Par::print(FILE *file) const {
	for(int i=0; i<NtabD; i++) { fprintf(file, " %12.9f", tabD[i]); }
	for(int i=0; i<NtabI; i++) { fprintf(file, " %6d", tabI[i]); }
}

void Par::print_verboten(FILE *file) const {
	fprintf(file, "PARAMETERS: n=%f, t=%f, J=%f, U=%f, alpha=%f, beta=%f, L=%d\n", n(), t(), J(), U(), alpha(), beta(), L());
}

void Par::print_label(FILE *file) const {
	fprintf(file, "n   t   tPrim   U   V   h   epsilonF  beta   tF   gFg    Vint    alpha    L");
}

// Var

void Var::print(FILE *file) const{
	for(int i=0; i<NtabD; i++) { fprintf(file, " %12.9f", tabD[i]); }
}

void Var::print_verboten(FILE *file) const{
	fprintf(file, "VARIABLES:");
	for(int i=0; i<NtabD; i++) { fprintf(file, " %22.15e", tabD[i]); }
	fprintf(file, "\n");
}

void Var::print_label(FILE *file) const {
	fprintf(file, "  lnF          lmF          nF             mu            mF           d       chiup        chidown      lchiup      lchidown");
}

// Res

void Res::set(gsl_multiroot_fsolver *s) {
	for(int i=0; i<n; i++) {
		(*V)[i] = gsl_vector_get (s->x, i);
	}
	this->glob = case_BF_sW::glob;
}

void Res::print(FILE *file) const { V->print(file); fprintf(file, "%22.15f %22.15f %12.9f %12.9f", glob.Energy, glob.nB, glob.qt, glob.qs); }
void Res::print_verboten(FILE *file) const { fprintf(file, "Result: "); V->print_verboten(file); fprintf(file, "%22.15f %22.15f %22.15f %22.15f\n", glob.Energy, glob.nB, glob.qt, glob.qs); }
void Res::print_label(FILE *file) const { V->print_label(file); fprintf(file, "             E                       nB             qs             qt             ");}

void example() {

    printf("Koniec\n");
}

void example_inf() {
	printf("Start BF infinity loop\n");

    Par initPar_mw(N_PAR);
    initPar_mw.set_n(0.8);
    initPar_mw.set_t(-1.);
    initPar_mw.set_J(1./3.);
    initPar_mw.set_U(14.);
    initPar_mw.set_alpha(0.5);
    initPar_mw.set_beta(500.);
    initPar_mw.set_L(256);
    initPar_mw.print_verboten();

	// I'm choosing start point
	Var initV(N_VAR);

	//FILE * pFile;
	//pFile = fopen("test2.txt", "w");

	// INFINITE
	Res res(N_VAR);
	RootGSL root(initPar_mw, function, res);
	printf("sdelta      "); res.print_label(); printf("\n"); fflush(stdout);
	while(1) {
	    initV.set_chiS(2.*frand()-1.);
	    initV.set_DeltaS(2.*frand()-1.);
	    initV.set_chiB(2.*frand()-1.);
	    initV.set_DeltaB(2.*frand()-1.);
	    initV.set_mu(8.*frand()-4.);
	    initV.set_lBF(8.*frand()-4.);

		if(root.Solve(initV) == GSL_SUCCESS) {
			printf("       "); root.print();printf("\n");  fflush(stdout);

		}
	}
	//fclose (pFile);
}

void run_1D() {


	return;
}


} // end of namespace case_BF_sW
